home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / statistics / lag1_source.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-05  |  1.7 KB  |  52 lines
  1. /* statistics/lag1_source.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Jim Davies, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20.  
  21. double 
  22. FUNCTION(gsl_stats,lag1_autocorrelation) (const BASE data[], const size_t stride, const size_t n)
  23. {
  24.   const double mean = FUNCTION(gsl_stats,mean) (data, stride, n);
  25.   return FUNCTION(gsl_stats,lag1_autocorrelation_m)(data, stride, n, mean);
  26. }
  27.  
  28. double 
  29. FUNCTION(gsl_stats,lag1_autocorrelation_m) (const BASE data[], const size_t stride, const size_t size, const double mean)
  30. {
  31.   /* Compute the lag-1 autocorrelation of a dataset using the
  32.      recurrence relation */
  33.  
  34.   size_t i;
  35.  
  36.   long double r1 ;
  37.   long double q = 0 ;
  38.   long double v = (data[0 * stride] - mean) * (data[0 * stride] - mean) ;
  39.  
  40.   for (i = 1; i < size ; i++)
  41.     {
  42.       const long double delta0 = (data[(i-1) * stride] - mean);
  43.       const long double delta1 = (data[i * stride] - mean);
  44.       q += (delta0 * delta1 - q)/(i + 1);
  45.       v += (delta1 * delta1 - v)/(i + 1);
  46.     }
  47.  
  48.   r1 = q / v ;
  49.  
  50.   return r1;
  51. }
  52.